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Figure 1: A “knot” model with a stripe texture (left) projected onto Chuang et aids function space (middle) and our 
connectivity-aware function space (right). Due to the coupling of function values, Chuang et al. ’s approach fails to repro¬ 
duce the correct texture when points are close in Euclidean space but are geodesically distant. By designing our function space 
to be aware of local connectivity, we can fit the geodesically distant patches independently, resulting in a better reproduction of 
the original texture. 


Abstract 

Recent work on octree-based finite-element systems has developed a multigrid solver for Poisson equations on 
meshes. While the idea of defining a regularly indexed function space has been successfully used in a number of 
applications, it has also been noted that the richness of the function space is limited because the function values 
can be coupled across locally disconnected regions. In this work, we show how to enrich the function space by 
introducing functions that resolve the coupling while still preserving the nesting hierarchy that supports multigrid. 
A spectral analysis reveals the superior quality of the resulting Laplace-Beltrami operator and applications to 
surface flow demonstrate that our new solver more efficiently converges to the correct solution. 

Categories and Subject Descriptors (according to ACM CCS): 1.3.3 [Computer Graphics]: Picture/Image 
Generation—Line and curve generation 


1. Introduction 

Solving the Poisson equation on meshes is essential in nu¬ 
merous geometry-processing applications. The task is most 


often formulated in terms of finite elements and two chal¬ 
lenges commonly arise: discretizing the space of functions 
on meshes and solving the resulting system of equations. 
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For discretizing the system, the most ubiquitous approach 
is to use tent-functions centered at mesh vertices [Dzi88], 
resulting in the classical cotangent-weight formula for tri¬ 
angle meshes [PP93]. For solving the sytem, black-box 
solvers like Algebraic Multigrid [RS87, HYOO] (iterative) 
and CHOLMOD [DH99,Davll] (direct) have been popular. 

Recently, Chuang et al. proposed an alternate framework 
that addresses both aspects simultaneously [CLB*09]. The 
idea is to define a function space in 3D and then restrict the 
3D functions to the mesh. There are several advantages in 
doing so: (1) the independence of the function space from 
the mesh helps define a Laplace-Baltrami operator that does 
not depend on the tessellation; (2) the nesting structure of 
the function space supports an efficient multigrid solver; and 
(3) the regularity of the function space can be leveraged in 
parallelizaing the solver. 

Though the approach has been successfully applied in a 
number of applications [CKlla, CKllb], it has been noted 
to produce artifacts because Euclidean distances are used to 
infer geodesic proximity. In particular, function values on 
locally disconnected components tend to be coupled when 
they are close in 3D, diminishing the richness of the function 
space (Figure 1). This also reduces the effectiveness of the 
multigrid solver, as disconnected regions are more likely to 
support the same basis function at coarser resolutions. 

The goal of this work is to address this coupling issue 
without sacrificing the regularity and nesting structure of the 
function space. The key idea is to make the function space 
connectivity-aware. The function space is enriched by split¬ 
ting existing functions such that the support of each new 
function is connected. 

The paper is organized as follows. After a brief literature 
survey in Section 2, we review Chuang et al .’s approach in 
Section 3, setting up a context that facilitates the presentation 
of our approach in Section 4. In section 5, we conduct a spec¬ 
tral analysis revealing the improved quality of our Laplace- 
Beltrami operator, we show the superior convergence of the 
resulting multigrid sovler, and we demonstrate the compet¬ 
itiveness of our solver with the state-of-the-art CHOLMOD 
solver [Davl 1] in a surface flow application. Finally, we con¬ 
clude by summarizing our work and discussing directions 
for future research in Section 6. 

2. Related Work 

Solving a Poisson-like system is a fundamental step in nu¬ 
merous geometry-processing applications, including mesh 
editing [SCOL*04], mesh deformation [SA07], surface re¬ 
construction [KBH06], surface fairing [Tau95, DMSB99], 
geometry sharpening/smoothing [CKllb], and surface pa¬ 
rameterization [FH05]. 

Defining the discrete Laplacian operator on meshes has 
attracted a great deal of research in the computer graph¬ 
ics community. Graph-based, combintorial operators have 


advantages of simplicity and efficiency [Tau95, Zha04]. 
Geometry-driven operators taking angles and areas into ac¬ 
count give rise the ubiquitous cotangent Laplacian [PP93, 
Flo03]. Recently, an operator based on intrinsic Delaunay 
triangulation has been proposed that addresses the prob¬ 
lem of non-convex weighting [FSBS06,BS07]. We refer the 
reader to [WMKG07] for a study of the tradeoffs between 
different operators. 

Solving the resulting system of equations is also of essen¬ 
tial importance. Poisson-like systems are usaully sparse and 
direct sparse solvers have been a popular choice when mem¬ 
ory is not a concern [DH99,Dav04,Li05]. In particular, when 
the system is symmetric and positive-definite, Cholesky fac¬ 
torization is often favored due to its numerical stability and 
efficiency [GVL96] . Furthermore, when the system is fixed 
but needs to be solved repeatedly for a constraint that varies 
over the course of an application, the factorization cost can 
be amortized, making these sparse factorization techniques 
an appealing option [SCOL*04]. 

Multigrid methods have also been used [BHMOO]. Their 
low memory usage and ability to dampen low-frequency 
errors efficiently make them preferable to direct solvers 
when the problem size is exceedingly large and/or the 
exact solution is not necessary. To define the hierarchi¬ 
cal structure on unstructured domains, “black-box” alge¬ 
braic multigrid methods [RS87] that rely solely on alge¬ 
braic information encoded in matrices have been studied 
[CFH*00,BCF*01,CFH*03]. Additionally, more geometry- 
driven approaches that rely on the design of mesh hierarchies 
have appeared in different contexts of geometry-procesing 
[KCVS98, CDROO, SK01, RL03, AKS05, SYBF06, SBZ09]. 

3. Review 

In this section, we briefly review the FEM formulation for 
solving the Poisson equations on meshes and Chuang et al.’s 
choice of function spaces [CLB*09, CKlla, CKllb]. 

3.1. Finite Elements on Meshes 

Given a Riemannian 2-manifold A4 and a continuous func¬ 
tion / : M R, solving the Poisson equations amounts to 
finding another function u on M whose Laplacian is /, i.e., 

A M U = f (1) 

with denoting the Laplace-Beltrami operator, the gener¬ 
alization of the Laplace operator to M. 

Instead of looking for the exact solution within the 
infinite-dimensional space of functions on Ad, a finite- 
elements system seeks an approximate solution within a 
finite-dimensional subspace B : A4 R spanned by the 
functions {bi(p), ■ • ■ ,b n (p)}. This is done by projecting 
both sides of Equation 1 onto B and solving an nxn lin- 
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ear system Lu = f where 


Lij = J ^bi(p)-bj(p)dp= -J('Vb i (p),Vb j (p))dp 

M M 

( 2 ) 

f i■ = J f{p)-bi{p)dp 

M 


The second equality in the first line uses the weak formula¬ 
tion derived from Stokes’ Theorem. It is correct when AT is 
closed or when either the test functions or the normal deriva¬ 
tives vanish at the boundary. The best solution within B is 
then u{p) = E'Lj Uibi(p). 


3.2. 3D B-Splines on Meshes 

To set up B , Chuang et al. start from a space of functions on 
R 3 (independent of Al) that are piecewise trilinear within a 
regular voxel grid Q N of resolution N xN x N (whose vox¬ 
els/corners we will denote by V(Q N )/K(Q N )). t 

Definition 1 Denote by B N the space of continuous func¬ 
tions in R 3 , with / G B N if and only if for each voxel 
v G V{Q N ) there exists a trilinear function f v : R 3 — R s.t. 
f(p) = Mp ) Vp G v. 

Property 1 (Nesting) B N C B 2N 

Proof When / G B N , Vv G V(G N ) 3 f v : R 3 R that is 
trilinear and f(p) = fv(p) Vp G v. Since Vv' G V ( Q 2N ) 3 v G 
V(Q n ) s.t. v' C v, it follows that f(p) = fv(p) Vp G v', and 
hence / G B 2N . □ 

Next, the function space B ^ on M is obtained by embed¬ 
ding AT within the voxel grid and considering the restriction 
of functions in B N to AT. 

Definition 2 Denote by B ^ the space of continuous func¬ 
tions on M, with / G B if and only if for each voxel 
v G V(G n ) there exists a trilinear function f v : R 3 — R s.t. 

f(p) = Mp) VpevHM. 

Property 2 (Nesting) B^ C B 2 ^ 

Proof When / G B N ^ Vv G V(Q N ) 3 f v : R 3 R that 
is trilinear and f(p) = f v (p) Vp G v fl AT. Since Vv' G 
V(Q 2N ) 3 v G V(Q n ) s.t. v' C v, it follows that f(p) = 
fv(p) Vp Gv' n AT, and hence / G B 2 ^. □ 

Using the nesting structure, Chuang et al. implement a multi¬ 
grid solver. Specifically, they use the trilinear B-splines {bf\ 
centered at grid comers k G K(Q N ) as test functions, and 
then define the prolongation operator by 

h(p)= I vt, N (k,k')-b k ,(p) (3) 

k'eK{Q 2N ) 


where V^ N is the standard prolongation stencil for 3D B- 
spline refinement obtained by taking the tensor-product of 
ID prolongation stencils: \ (1331). 


4. Approach 

A limitation of Chuang et aV s approach is that the notion 
of continuity is characterized with respect to the Euclidean 
distance rather than the geodesic one. One manifestation 
of this is that two points adjacent in 3D will necessarily 
have similar function values, even if 
they are geodesically distant. The in¬ 
set shows a 2D example, where each 
of the four B-splines supported on the 
cell has similar values on p and q. As 
a result, the values of any function in 
B 1 ^ will be coupled at p and q. 

4.1. A New Function Space 

The goal of this work is to address the coupling problem by 
enriching B ^. In particular, if there exist multiple connected 
components within a voxel, we will consider assigning a sep¬ 
arate basis function to each of them. 

Definition 3 We define B ^ as the space of continuous func¬ 
tions on AT with / G B ^ if and only if for each voxel v G 
V(G N ) and each connected component c C v D AT, there ex¬ 
ists a trilinear function f c : R 3 —)• R s.t. f(p) = fc(p)Vp £ c. 

Observation If c is a connected component of v' fl AT for 
some voxel vGf ( Q 2N ) and c is a connected component of 
v n AT for some voxel v G V ( Q N ), then either c fl c = 0 or 
c C c. Additionally, for each c C v' G V(Q 2N ), there must 
exist acCvGV (G N ) s.t. v' Cv and c C c. 

Claim (Nesting) B^ C B 2 ^ 

Proof When / G By^ , Vv G V{Q N ) and each connected com¬ 
ponent cCvnAT,3/ c :R 3 —>>R that is trilinear and f(p) = 
f c (p) V/?Gc. Since Vv' G V(G 2N ) and each connected com¬ 
ponent c CvH AT, 3c C v fl AT for some v G V(Q N ) s.t. 
c C c, it follows that f(p) = f c (p) Vp G c', and hence 
feB™. □ 

Note that, as^with the work of Chuang et al ., the nesting 
structure of B supports the definition of a multigrid solver. 

The advantage of this formulation is highlighted in Fig¬ 
ure 1 , where the texture on the self-intersecting mesh on the 
left is projected onto B ^ and B%_. In this example, there 
does not exist a continuous 3D function whose restriction to 
the mesh can closely represent the texture. 



4.2. Implementation: Finding The Basis 

t Though in this work we focus on the piecewise trilinear space of Chuang et al. use first-order, tensor-product B-splines po- 

functions, higher order spaces can also be set up, as with [CLB*09]. sitioned at corners of a regular grid and then discard those 
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Eigenvalue Eigenvalue 




Figure 3: Stability of the spectrum of the estimated Laplace-Beltrami operator We compute the spectra of Chuang et al.’s 
operator (left) and our operator (right) at various grid resolutions for the Pulley model (Figure 4, bottom left). As the resolution 
increases, our spectra more quickly converge to the ground-truth. 



Figure 2: Adaptive splitting of test functions. In [CLB*09], 
the test functions are chosen independent of the mesh, result¬ 
ing in disconnected components in the support (left). In con¬ 
trast, our approach refers to mesh connectivity and assigns a 
separate test function to each component (middle and right). 


test functions whose support does not overlap the surface, to 
obtain the test functions spanning . 


In order to obtain our test functions, we start by inspect¬ 
ing the supports of Chuang et aids functions { b k }. At each 
corner k E K(Q N ), we generate a separate test function for 
each connected component of supp(^) f] M (that is, the in¬ 
tersection of M with the 2x2x2 voxels around k). This is 
illustrated in Figure 2 using a 2D example. We will denote 
by 4 the number of connected components in the support of 
b k , and denote by C l k C M the j-th connected component. 
Then, each of our test functions b l k is duplicated from b k but 
is supported only on C[: 


b’k(p) = 


h(p) if P £ C[ 

0 otherwise. 


4.3. Implementation: Multigrid 

We obtain our multigrid hierarchy by setting up the test func¬ 
tions (as described in Section 4.2) using grids of successively 
finer resolutions. As in Chuang et al. ’s work, every test func¬ 
tion b\ defined by grid Q N can be expressed as the linear 


combination of the test functions {b l k ,} defined by grid Q 2N . 
Specifically, the prolongation operator is defined by 

bi(p)= £ t d x(clclfvl N {k,k')-liAp) (4) 

k'eK(g 2 N )i '=i 

where 

y ( CC ') = l 1 ifcnC '^0 

M 5 ; \ 0 otherwise 

indicates if the component c is contained in c. 


5. Results 

In this section, we describe several experiments for evaluat¬ 
ing our approach. We start with a spectral analysis that re¬ 
veals the robustness of our Laplace-Beltrami operator. Next, 
we examine how the rate of solver convergence is improved 
by our connectivity-aware space of functions. Finally, we ap¬ 
ply our approach to a surface flow application and compare 
its performance with the state-of-the-art CHOLMOD solver. 


5.1. Spectral Analysis 

Eigendecomposition of the Laplace-Beltrami operator has 
been widely used for analyzing and processing signal on 
meshes [Tau95, VL08]. 

To evaluate the quality of the estimated Laplace-Beltrami 
operator (Equation 2), Chuang et al. compare the spectrum 
of their operator with the spectrum of the cotangent-weight 
operator. In particular, this is done by solving the generalized 
eigenvalue problem: 

Lx = XMx 

with x a generalized eigenvector and M the mass matrix 

Mij= / bi(p)-bj(p)dp 
J M 
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Rotated Models 





Figure 4: Rotation invariance of the estimated Laplace-Beltrami operator. We compute the spectra of Chuang et al.T opera¬ 
tor (top and middle) and our operator (bottom) for different rotations of the Pulley model. The zoom-ins accentuate the superior 
stability of our operator (right). 


They show that the idea of restricting regular 3D functions 
generally helps define a robust Laplace-Belrami opearator 
whose spectrum is more stable even at a lower dimension 
(we refer the reader to [CLB*09] for details). However, as 


they also point out, the approach has difficulty in regions 
with “narrow cross-sections” - precisely where the coupling 
of values occurs. 

We rerun the experiment on the same model with which 
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Figure 5: Color fitting of a 3D checker-board texture on the model consisting of equidistantly-spaced 6x6x6 unit-cubes 
(left). The screened-Poisson equation with the screening weight a = 0.01 is solved using a grid of depth 5 (i.e., consisting of 
2 5 x 2 5 x 2 5 voxels). The coefficients of the initial guess are generated randomly with values between 0 and 1. The residuals 
(normalized by the norm of the initial guess) after one W-cycle are ploted as a function of the minimum depth of the multigrid 
solver (right), ranging from 0 (i.e., full W-cycle relaxation) to 5 (i.e., Gauss-Seidel relaxation at the finest level only). The texture 
is reconstructed using Chuang et al. ’s approach (middle left) and our connectivity-aware approach (middle right). 


they had trouble (Figure 4, left). We compute the spectrum 
at increasing grid resolutions (Figure 3). As demonstrated in 
the plot, the spectra of our operator quickly converge to the 
ground truth (computed using the cotangent-weight operator 
defined over a dense tessellation of the mesh). 

We also verify that the rotation-invariance property of the 
Laplacian is preserved (Figure 4). We apply different rota¬ 
tions to the model before computing the spectrum. Note that 
even when defining Chuang et al.’s operator using a higher 
resolution grid, our operator still reproduces the spectrum 
over the different rotations more consistently. 


5.2. Solver Convergence 

As shown in Figure 1, sometimes the connectivity-unaware 
function space is just not rich enough to express the desired 
solution. In this section, we examine how the coupling also 
affects the convergence of the multgrid solver and show that 
our connectivity-aware approach helps resolve the issue. 

In order to make the solver residuals from the different 
spaces comparable, we use an input model (Figure 5, left) 
which defines the same function spaces B 1 ^ and B at suf¬ 
ficiently fine resolution N. This means that the splitting oper¬ 
ation of our approach is not necessary at the finest resolution 
(though it is still necessary at coarser ones). We then solve a 
screened-Poisson equation on the mesh to reconstruct a color 
function /: M —» M 3 . That is, given the function /, we look 
for the coefficients u satisfying: 

(L + a • M) u = (f + a • s) 


with constraints 

f;= f (Vf(p),Vbi(p))dp 

JM 

Si= [ f(p)'bi(p)dp 
J M 

and screening weight a = 0.01. Note that the screening here 
serves to softly constrain the constant term on each con¬ 
nected component. We perform one standard W-cycle relax¬ 
ation with ten iterations of Gauss-Seidel smoothing at each 
level. All coefficeints are intialized with random values be¬ 
tween zero and one. 

From the results of the reconstruction (Figure 5, middle), 
we observe that the connectivity-unaware solver does not 
produce a satisfactory solution, despite the lack of coupling 
at the finest resolution. We believe this is because the cou¬ 
pling at the coarser spaces prevents multigrid from generat¬ 
ing a meaningful correction term. 

The plots on the right confirm this. Here we show the 
magnitude of the residuals as the minimum depth of the 
multigrid solver is increased. Looking at this plot, we see 
that increasing the minimum depth does not affect the con¬ 
vergence of the connectivity-unaware solver, indicating that 
the coarser correction terms do not improve the solution. On 
the other hand, the convergence of our connectivity-aware 
solver quickly deteriorates when we increase the minimum 
depth, suggesting the correct multigrid behavior. 

5.3. Surface Flow 

As pointed out by Chuang et al. , although the grid-based ap¬ 
proach helps define an effective multigrid, the preprocessing 
time is significant compared to the solver time. This makes 
the approach particularly suitable for applications where an 
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Figure 6: Conformalized Mean-Curvature Flow applied to the Armadillo Man. From left to right, we show the Oth, 1st, 3rd, 
5th, 10th and 30th steps of the flow. The flow is numerically stable and conformally evolves the mesh to a sphere (right). 


evolving linear system needs to be solved repeatedly. We 
evaluate our method in this context. 

We apply our approach to support conformalized mean- 
curvature flow (cMCF) [KSBC12]. The flow, recently pro¬ 
posed by Kazhdan et al. , has been shown to converge to 
a (conformal) spherical parameterization when acting on 
genus-zero surfaces. Importantly, as opposed to traditional 
mean-curvature flow, cMCF does not develop singularities 
and is numerically stable, making it well-suited for studying 
the long-term behavior of our solver. 

At each time t , we solve a semi-implicit system as de¬ 
scribed in [DMSB99]: 

5 A 

Mt F-L 0 J u t+5 = M t ut 

with 5 the temporal stepsize, M t the mass matrix of the sur¬ 
face at time t , and Lq the stiffness matrix of the original sur¬ 
face. Figure 6 shows an example of the flow. 


Performance When performing surface flow, it is often 
necessary to consider the trade-offs between the computa¬ 
tional cost and the solution accuracy. Factors affecting the 
computational cost include the temporal stepsize (as taking 
smaller timesteps increases the temporal resolution but leads 
to longer running time) and the solver time per step (as us¬ 
ing more accurate solvers increases the accuracy within each 
timestep but leads to longer running time). 

We first simulate the ground truth cMCF on a high reso¬ 
lution brain model consisting of 1.4 millions vertices (Fig¬ 
ure 8, left). The evolution time is targeted at t = 50 and we 
take a tiny stepsize 8 = 0.01 to flow the surface toward the 
target. At each time step, we use the cotangent-weight op¬ 
erator to define the system, which is then solved precisely 
by CHOLMOD. The simulation takes more than 10 hours 
to generate all the evolved surfaces on a PC with an \1- 
2820QM processor. Having computed these, we later gen¬ 
erate the ground truth at any t = T by linearly interpolating 
the surfaces at t = |_|J5 and t = [|]8. 


Next, we run the flow using the following configurations 
of systems and solvers: * 

• Cotangent-weight System solved by CHOLMOD 

• Chuang et al.’s system solved by CHOLMOD 

• Chuang et al.’s system solved by Multigrid 

• Our system solved by CHOLMOD 

• Our system solved by Multigrid 

For the grid-based systems (all but the first), we follow 
Chuang et al.’s formulation of test function tracking, in or¬ 
der to avoid having to set up the multigrid hierarchy repeat¬ 
edly as the embedding changes [CKlla]. We choose a grid 
resolution of 2 8 x 2 8 x 2 8 so that the resulting dimension is 
roughly the same as that of the cotangent-weight system. § 

We require each configuration to complete the flow within 
one hundred seconds. As the evolution time is fixed at t = 50, 
the temporal stepsize 8 taken by each configuration depends 
on how quickly each spacial system is solved, as summa¬ 
rized in Table 1 . 

Finally, we compare the evolved surfaces obtained from 
each configuration to the ground truth. Results are shown 

in Figure 7, where we plot the RMS error (y^ jjv? — vf || 2 
with v e and v g the evovled and ground truth vertex posi¬ 
tions) as a function of evolution time. Here we make two 
observations. First, although the direct solver is capable of 
computing exact solutions, the expensive cost prevents us 
from taking small timesteps and eventually leads to larger 
errors. Second, the connectivity-unaware and connectivity - 
aware multigrid systems perform equally well in the be- 
gining, but then the connectivity-unaware one deteriorates 
quickly with time. In Figure 8, we examine the two surfaces 
at the end of the flow. The zoom-ins highlight the problem of 
the connectivity-unaware approach, where the values of the 


t Our CHOLMOD is compiled by the Intel Math Kernel Li¬ 
brary [MKL12] to ensure the best performance. In addition, the 
symbolic factorization is only performed once in the preprocessing. 
§ When the dimensions match, solving the grid-based systems is 
more costly than solving the cotangent-weight system, as the grid- 
based systems are less sparse (the average number of entries per row 
is 18, compared to 7 for the cotangent-weight Laplacian). 
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Configuration 
(System + Solver) 

Dimension 

Non-Zero Entries 

Time/Step 

Steps 

Step size (5) 

Cotangent + CHOLMOD 

1.38 x 10 6 

9.67 x 10 6 

3.78 

26 

1.92 

+ CHOLMOD 



11.08 

9 

5.56 

[CLB* 09] wl . 

+ Multigrid 

1.21 x 10 6 

2.22 x 10 7 

1.47 

68 

0.74 

^ + CHOLMOD 



9.26 

10 

5.00 

Ours wl . ., 

+ Multigrid 

1.27 x 10 6 

2.28 x 10 7 

1.52 

66 

0.76 


Table 1: Statistics for the different configurations, giving the system dimension, the number of non-zero entries, the average 
time spent for each time step (including the time for updating the system/solver and the time for solving the system), the total 
number of steps, and the temporal step size. 


coordinate function are coupled across the two hemispheres 
of the brain and cannot flow independently. 


Discussion 

Using our formulation, we resolve the problem of value cou¬ 
pling across disconnected components. However, we cannot 
decouple function values for points on the same component. 
The problem becomes manifest when the grid resolution is 
low and the surface has high curvature. In Figure 10, we vi¬ 
sualize the support of one basis function of our system de¬ 
fined using a lower resolution grid (depth 7). The support 
contains two flat regions that meet near a corner. As a result, 
running cMCF at this low resolution results in pronounced 
errors in these regions. 

Another issue with grid-based systems, is the possible 
presence of linear dependency. That is, when there exists 
an entirely planar component aligned with one of the axes 
(i.e. not in general position ), the linear system becomes sin¬ 
gular. In pratice this is not a problem, because our multigrid 
solver uses a Gauss-Seidel smoother, which is guaranteed 
to converge for symmetric and positive definite systems and 
tends to behave well even when the system is only positive 
semi-definite. This can become a problem when using a di¬ 
rect solver like CHOLMOD. However, this issue can be re¬ 
solved, either by applying a slight rotation to the model, or 
by adding a diagonal term, 8 • Id, with 8 a small constant, to 
the system. 

Finally, one can interpret our connectivity-aware adap¬ 
tation of the FEM construction as a refinement of a 
point-set topology on the surface. Specifically, if one de¬ 
fines a simple topology on R consisting of the open sets 
X(R) = {0,R — {0},R}, then any choice of basis functions 
{b\ ,..., b^} defines a topology X(M) on the mesh, gener¬ 
ated by the unions and intersections of the sets bf 1 (U), with 
U E X(R). By splitting basis functions based on connectiv¬ 
ity, we obtain a finer set set of generators for a topology 
on M, reflecting the better localization of the connectivity- 
aware function space. 


6. Conclusion 

In this work, we present an important extension to grid-based 
systems for solving Poisson equations on surfaces. Our for¬ 
mulation addresses the coupling issue arising in earlier work 
by incorporating local connectivity information in defining 
an FEM discretization. The resulting function space has sev¬ 
eral advantages: It defines a robust Laplace-Beltrami oper¬ 
ator (5.1). It improves the multigrid behavior (5.2). And, it 
provides an effective system for performing surface evolu¬ 
tion (5.3). 

In the future, we would like to explore extension of our 
approach to feature-adaptive refinement of test functions (in 
particular, near the high-curvature regions where the cou¬ 
pling problem still occurs). We would also like to apply 
our approach to support other applications that can benefit 
from a real-time iterative solver, such as cloth simulation and 
more general surface flows. 
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Figure 7: Error comparison of the different approaches performing cMCF on a brain model consisting of 1.4 million vertices. 
The RMS error is plotted as a function of evolution time. The ground truth is simulated using CHOLMOD to solve the cotangent- 
weight system taking a tiny stepsize 5 = 0.01. The computational budget is fixed at one hundred seconds, so that the number of 
steps (visualized by the tick marks) is determined by the efficiency of the solver. 



Figure 8: The brain undergoing cMCF with the connectivity-unaware ([CLB* 09]) and connectivity-aware (ours) function 
spaces used to define the system for the input mesh (left). The two systems have about the same dimension and are both solved 
using the multigrid. Here we show the evovled surfaces at t = 10, t = 25, and t = 50 (middle). The zoom-ins highlight the benefit 
of using the connectivity-aware system, which is able to decouple the function values at points that are close in Euclidean space, 
allowing them to flow independently (right). 
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Figure 9: Error visualization of the cMCF evolution by the different configurations at time t = 50. We draw on the mesh 
the per-vertex squared distance to the ground truth in red. Note that for the connectivity-unaware systems, errors accumulate 
quickly around those “ pinched ” regions. 



Figure 10: Value coupling within a connected component near a high curvature region of the brain (left). When we define our 
system using a lower resolution grid, there exists a basis function supported on two parallel patches. Rendering the support 
from different perspectives, we observe that the function cannot be split because the two partches are connected at the corner 
(middle). As a result, performing cMCF using this lower resolution system yields high errors on the evolved surface (left, drawn 
as in Figure 9). 
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